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HIGHLIGHTS 


• Our Models resulted in very high out-of-sample R 2 ( R 2 = 0.88). 

• Our model performed well both spatially and temporally (R 2 = 0.87, R 2 = 0.87). 

• Our results revealed very little bias (Slope of predictions versus withheld observations = 0.99). 

• Importantly, these R 2 are for daily, rather than monthly or yearly, values. 
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The use of satellite-based aerosol optical depth (AOD) to estimate fine particulate matter (PM 2 . 5 ) for 
epidemiology studies has increased substantially over the past few years. These recent studies often 
report moderate predictive power, which can generate downward bias in effect estimates. In addition, 
AOD measurements have only moderate spatial resolution, and have substantial missing data. We make 
use of recent advances in MODIS satellite data processing algorithms (Multi-Angle Implementation of 
Atmospheric Correction (MAIAC), which allow us to use 1 km (versus currently available 10 km) reso- 
lution AOD data. We developed and cross validated models to predict daily PM 2.5 at a 1 x 1 km resolution 
across the northeastern USA (New England, New York and New Jersey) for the years 2003-2011, allowing 
us to better differentiate daily and long term exposure between urban, suburban, and rural areas. 
Additionally, we developed an approach that allows us to generate daily high-resolution 200 m localized 
predictions representing deviations from the area 1 x 1 km grid predictions. We used mixed models 
regressing PM 2.5 measurements against day-specific random intercepts, and fixed and random AOD and 
temperature slopes. We then use generalized additive mixed models with spatial smoothing to generate 
grid cell predictions when AOD was missing. Finally, to get 200 m localized predictions, we regressed the 
residuals from the final model for each monitor against the local spatial and temporal variables at each 
monitoring site. Our model performance was excellent (mean out-of-sample R 2 = 0.88). The spatial and 
temporal components of the out-of-sample results also presented very good fits to the withheld data (R 2 
= 0.87, R 2 = 0.87). In addition, our results revealed very little bias in the predicted concentrations (Slope 
of predictions versus withheld observations = 0.99). Our daily model results show high predictive ac- 
curacy at high spatial resolutions and will be useful in reconstructing exposure histories for epidemio- 
logical studies across this region. 

© 2014 Elsevier Ltd. All rights reserved. 
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1. Introduction 

The use of satellite-based aerosol optical depth (AOD) to esti- 
mate exposure to fine particulate matter (PM 2 . 5 -particulate matter 
<2.5 pm in aerodynamic diameter) concentrations for epidemio- 
logic studies has increased substantially over the past few years 
(Chang et al., 2013; A.A. Chudnovsky et al., 2013; Chudnovsky et al., 
2012 ; Kim et al., 2013 ; Kloog et al., 2012b, 2011 ; Lee et al., 2011 ; Lin 
et al., 2013; Nordio et al., 2013). Both short term (hours, days) and 
chronic (months, years) exposure to PM 2.5 has been extensively 
associated with detrimental human health effects (Halonen et al., 
2008; Kloog et al., 2012a; Peacock et al., 2011; Zanobetti and 
Schwartz, 2009). Fine PM exposures were found to be associated 
with respiratory and cardiovascular morbidity, mortality from 
cardiovascular and respiratory diseases and an increase in hospital 
admissions (Dominici et al., 2006; Kloog et al., 2012a; Schwartz, 
1996). 

Traditionally, estimation of PM2.5 exposures has been based on 
using measurements of a central ground monitor, and assigning 
these measurement to populations within a specified distance of 
the monitor (Laden et al., 2006; Samet et al., 2000). This approach 
introduces exposure error (also known as exposure misclassifica- 
tion), and likely biases the effect estimates downward due to spatial 
misalignment (Zeger et al., 2000). In recent years, to avoid this 
exposure misclassification, many studies have used regression 
models based on geographic covariates (termed “land use” re- 
gressions — LUR) to expand in situ measurements of PM2.5 con- 
centrations to large areas. LUR is essentially an interpolation 
technique that employs the PM2.5 concentrations as the dependent 
variable, with proximate land use, traffic and physical environ- 
mental variables used as independent predictors (Beckerman et al., 
2013; Gryparis et al., 2009; Hoek et al., 2008; Vienneau et al., 2010). 
Since the geographic covariates are generally not time varying, the 
temporal resolution of the LU model predictions is limited and thus 
can only be employed to assess long-term exposures for chronic 
health effects studies. Finally, current exposure assessments are 
limited to urban areas due to the paucity of monitoring sites in rural 
locations. 

Satellite-based AOD provides physical daily measurements that 
can be used to estimate air quality and pollution due to its extensive 
spatial coverage and repeated observations of the earth surface and 
atmosphere. AOD is a measure of the extinction of electromagnetic 
radiation at a given wavelength due to the presence of aerosols in 
an atmospheric column. However, satellite-based AOD is a measure 
of light attenuation in the column which is affected by ambient 
conditions (such as vertical profile, chemical composition and hu- 
midity). In contrast, PM2.5 is a measure of dry particle mass near the 
surface and thus we do not expect them to be strongly correlated 
(Chudnovsky et al., 2012). 

In recent years many studies have used various statistical 
methods to establish quantitative relationships between AOD and 
PM2.5 (Chang et al., 2013; A. Chudnovsky et al., 2013; Cordero et al., 
2013; Gupta et al., 2013; Kim et al., 2013). Traditionally, the health 
exposure studies have used the standard MODIS (Moderate Reso- 
lution Imaging Spectroradiometer) AOD product of the “Dark 
Target” algorithm (Levy et al., 2007) which has a nadir resolution of 
10 x 10 km 2 . Lately, AOD at significantly higher spatial resolution 
(lxl km 2 ) has been offered by a new Multi-Angle Implementation 
of Atmospheric Correction (MAIAC) algorithm. Chang et al. (2013) 
used some novel ideas such as a statistical downscaling and data 
fusion techniques to predict PM2.5 concentrations at spatial point 
locations in the southeastern United States during the period 
2003—2005. Their model showed relatively high cross-validated 
predictions (R 2 = 0.78 and a root mean-squared error (RMSE) of 


3.61 |ig/m 3 ). Chudnovsky et al. (2014) used one year of observations 
from MODIS based Aqua at a 1 km spatial resolution to obtain daily 
PM 2.5 estimates for AOD retrieval days. Their model results indi- 
cated that high resolution MAIAC (Multi-Angle Implementation of 
Atmospheric Correction) data better explains the variations in 
PM 2.5 However, the developed model was restricted to retrieval 
days, and the out-of-sample predictive performance spatially was 
similar to that estimated by Chang et al. (spatial R 2 = 0.79). 

These recent studies, while generally showing better fits than 
previous models, still leave room for improvement in terms of 
reducing exposure error. In addition, they all lack detailed high 
resolution predictions across large space-time domains (especially 
in rural areas) which is critical for acute exposure epidemiologic 
studies. 

In this paper, we incorporate MAIAC-based AOD satellite data 
which allows us to present a much simplified model, while still 
improving significantly over currently available models. We devel- 
oped and validated models to predict daily PM 2.5 at a 1 x 1 km 
resolution across the northeastern USA (New England, New York 
and New Jersey) for the years 2003—2011, allowing us to better 
differentiate daily and long term exposure between urban, subur- 
ban, and rural areas. Additionally, we developed an approach that 
allows us to generate daily high-resolution 200 x 200 m localized 
predictions that are separate from the area lxl km grid predictions. 

2. Material and methods 

2.1. Study domain 

The spatial domain of our region included the Northeastern part 
of the USA (Fig. 1 ), and includes the states of Connecticut, Maine, 
Massachusetts, New Hampshire, New Jersey, New York, Rhode Is- 
land and Vermont. Many urban areas (notably New York and Bos- 
ton) are included as well as rural hill towns, large forested regions, 
water bodies, mountains and the Atlantic sea shoreline. The study 
area included 285,284 discrete lxl km satellite grid cells. 

2.2. AOD data 

One of the fundamental aerosol products from MODIS is spectral 
AOD (sometime referred to as Aerosol Optical Thickness or AOT). 
This is a global level 2 daily MODIS product. The details of the 
standard “Dark Target” (DT) MODIS algorithm over land providing 
10 x 10 km 2 resolution AOD have been broadly reported (Levy et al., 
2007; Remer et al., 2005). Recently a new processing algorithm 
(MAIAC) has been developed that provides a high 1 km resolution 
AOD from MODIS data (Chudnovsky et al., 2014; Lyapustin et al., 
2011a, 2011b). MAIAC processing begins with the gridding of 
MODIS LIB data to a fixed 1 km grid, and accumulation of up to 16 
days of measurements in the memory. Then it uses time series 
analysis and processing of groups of pixels to derive surface bidi- 
rectional reflectance distribution function (BRDF) and aerosol pa- 
rameters without assumptions typical of current MODIS 
operational processing algorithms. The spatio-temporal analysis 
also helps MAIAC's cloud mask augmenting traditional pixel-level 
cloud detection techniques. In this analysis we used MAIAC AOD 
based on collection 6 MODIS Aqua LIB data for the years 
2003-2011. 

2.3. Monitoring data 

Data for daily PM2.5 mass concentrations across the Northeast 
region (see Fig. 1 ) for the years 2003—2011 were obtained from the 
U.S. Environmental Protection Agency (EPA) Air Quality System 
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Fig. 1 . Map of the study area showing the location of EPA and IMPROVE monitoring sites across the Northeastern USA. 


(AQS) database as well as the IMPROVE (Interagency Monitoring of 
Protected Visual Environments) network. IMPROVE monitor sites 
are located in national parks and wilderness areas while EPA 
monitoring sites are located across the Mid-Atlantic including ur- 
ban areas such as New York City, Boston, etc. There were 161 
monitors with unique locations operating in the Northeast during 
the study period. The Mean PM 2.5 across the Northeast during the 
study period was 11.7 pg/m 3 with a standard deviation of 7.8 pg/m 3 
and an interquartile range (IQR) of 9.1 pg/m 3 . 

2.4. Spatial predictors of PM 2.5 

We used a hybrid model, containing both AOD, land use, and 
meteorological predictors. The following spatial predictors were 
used: population density, elevation, traffic density, percentages of 
land use according to 12 land use categories based on United States 
Geological Survey (USGS) National Land Cover Dataset (NLCD), 
point emissions and total area-source emissions (tons per year) for 
PM 2 . 5 , PM 10 , SO 2 , and NO*. 

Elevation : Elevation data were added through a satellite-based 
digital elevation model from the USGS National Elevation Dataset 
(NED) covering the United States at a spatial resolution of 1 arc sec 
(Maune, 2007). There are sharp elevation contrasts across such a large 


study area and thus we used elevation as a spatial predictor (generally 
higher elevations are associated with lower air temperatures). 

Traffic density : Road data were obtained through the U.S. Census 
2000 topologically integrated geographic encoding and referencing 
system (TIGER, 2006). We calculated the total Al road length (class 
1 roads that are hard surface highways including Interstate and U.S. 
numbered highways, primary State routes, and all controlled access 
highways) across the study area. The Al roads were intersected 
with the lxl km grid and the resulting attribute tables contained 
the density of all Al road segment lengths in each 1 km 2 grid cells. 

Population density: Population density data were obtained 
through the U.S. Census 2000 dataset (Census, 2000). We calculated 
the weight-averaged population for each lxl km grid cell based 
on the tracts intersecting these grid cells. 

Point emissions : Additional point emissions data for PM 2 . 5 , PM 10 , 
SO 2 , and NO* were obtained through the 2005 U.S. EPA National 
Emissions Inventory (NEI) facility emissions report (EPA, 2010). 
Locations reporting zero emissions within the appropriate grid cell 
were assigned a value of one-half of the minimum value among all 
monitoring locations. 

Area-source emissions: Area-source PM 2 . 5 , PM 10 , SO 2 , and NO* 
emissions data were obtained through the 2005 U.S. EPA-NEI tiered 
emissions reports (EPA, 2010), which provide estimates of total 
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area-source emissions by county and year. Intersecting source 
emission areas for each lxl km grid were calculated by weight 
averaging the values from each intersecting areas intersecting n 
each lxl km grid cell. 

Percentages of land use : We used the NLCD from 2001 (Homer 
et al., 2004), available as raster files with a 30 m spatial resolu- 
tion. We calculate the percentages of mixed forest, deciduous for- 
est, evergreen, crop, pasture, grass, shrub, water, high 
development, medium development, low development and open 
development areas in each lxl km grid across the study area. 

2.5. Temporal predictors of PM 2.5 

Meteorological data : All meteorological variables used in the 
analysis were obtained through the National Climatic Data Center 
(NCDC, 2010). Only continuous operating stations with daily data 
running from 2003 to 2011 were used (26 stations spread across the 
study area). Grid cells were matched to the closest weather station 
with available meteorological variables (24 h means). 

We used the following meteorological variables: air tempera- 
ture, wind speed, daily visibility, sea land pressure (SLP) and rela- 
tive humidity. All the variables represent 24 h averages except 
visibility, which is only computed during daylight hours. 

NDVI: We used the publicly available monthly MODIS NDVI 
(Normalized Difference Vegetation Index) product (MOD13A3) at 
1 km spatial resolution. The monthly resolution was chosen since 
NDVI values do not change considerably within a month except 
periods of spring green-up and fall senescence. 

PEL : We used publicly available daily data on the height of the 
planetary boundary layer (PBL) obtained from the NOAA Reanalysis 
Data (NOAA, 2010). The spatial scale of the data was 32 x 32 km. 
The height of the boundary layer may vary with wind speed (Oke, 
1987), influencing the concentration and vertical profile of pollut- 
ants. The boundary layer not only controls transport and location of 
pollutants and aerosols but also their concentrations would be 
different in variable boundary layer structures (Angevine et al., 
2013). 


geocoded data (RGD). When only SAGD health data is available, we 
use predictions at the lxl km grid cell level, whereas when RGD 
are available, we add an additional local daily estimation compo- 
nent at a very high resolution (200 x 200 m) to the grid-level 
predictions. Using higher resolution MAIAC data compared to our 
previous 10 x 10 km MODIS AOD data allowed us to simplify our 
approach in this way. 

All models were fit to data from each year (2003—2011) sepa- 
rately. To generate the daily lxl km PM2.5 predictions in each grid 
cell for the entire 2003—2011 period, we developed a prediction 
process using a series of models: The first model calibrates the AOD 
grid-level observations to the PM2.5 monitoring data collected 
within 1 km of an AOD value, while adjusting for the land use and 
meteorological variables. In the second model we predict daily 
PM2.5 concentrations in grid cells without monitors but with 
available AOD measurements using the model 1 fit. Then in the 
third model, to estimate daily PM2.5 in grid cells with no AOD on 
that day, we take advantage of the region-specific association be- 
tween grid-cell AOD and PM2.5 levels, and the association between 
PM2.5 level in a given grid with that in neighboring grid cells. 

To generate the daily 200 x 200 m PM2.5 local predictions for 
studies using RGD, we take the residuals of model 3 at each 
monitoring site (that is the ground level PM2.5 data minus the 
model 3 predictions) and regress them against very fine 
(200 x 200) monitor-specific spatial and temporal predictors of 
PM2.5 which included traffic density, population density, elevation, 
percent urban, distance to major roads (Al), distance to source 
emission points, PBL and visibility. This stage allows us to predict 
local PM throughout the study area on any given day. 

To accommodate the fact that daily AOD data missingness is not 
random, the first stage model incorporated inverse probability 
weighting (IPW) to potentially avoid bias in the regression coeffi- 
cient estimates and thus in the resulting predictions. This approach 
effectively up-weights dates and grid cells which are under- 
represented due to a large degree of missing data. To obtain the 
weights that account for the non-random missingness in AOD 
values, we fit the following logistic regression model for the 
probability (p) of observing an AOD value in cell i on day j: 


2.6. Statistical methods 

All modeling was done using the R statistical software version 
3.02 and SAS (Statistical Analysis System) version 9.3. 

As we have shown in previous studies (Chudnovsky et al., 2014; 
Kloog et al., 2012b, 2011) there is a varying spatial relationship 
between AOD and PM2.5 on a single day due to differences in factors 
including particle composition, PBL, relative humidity, vertical 
profiles. Thus, we use a mixed effects model incorporating spatial 
and temporal predictors and day-specific random-effects to take 
into account these temporal variations in the PM2.5— AOD rela- 
tionship. Since the exposure predictions generated by our model 
are primarily intended for use in epidemiologic health effect 
studies, we generate exposure predictions at two different spatial 
scales at which health data are typically collected: small area (zip 
code, census tract, etc) geocoded data (SAGD) and address-specific 


In 


P 

1 ~P 


= b 0 + ^Elevation* + b 2 SLPy + b 3 Temperaturey 
+ Month 


IPW= i (1) 

where (p) is the probability for availability of AOD in each day in 
each grid cell in each year. There were no observations which had a 
disproportionate influence in the yearly models. 

Finally, since the daily PM— AOD calibration factors can vary 
spatially, the study area is divided into 7 regions. The day-specific 
intercept, AOD, and temperature random effects in the model are 
nested within regions of the study. Specifically, the first model can 
be written as 


5 28 2 

PMy = a + Uj + gj (reg) + bi + Vj + hj(reg) AODy + b 3 + l<j Temperaturey + 9 1 m^l mi + Q 2 m^ 2 mj + 93m^3my + ij 

m= 1 m= 1 m= 1 

tljVjkj N[(000);s] 

Sj( reg) reg) [(00) ; s REg] 


(2) 
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where PM,) is the measured PM2.5 concentration at a spatial site i on 
dayj; a and uj are the fixed and random (day-specific) intercepts, 
respectively, AODy is the AOD value in the grid cell corresponding to 
site i on day j; bi and vj are the fixed and day-specific random 
slopes, respectively. Temperature,;/ is the temperature value in the 
grid cell corresponding to site i on day j; and, b2 and kj are the fixed 
and random slopes for temperature. X lm ,- is the value of the mth 
spatial predictor at site i, X 2rry - is the value of the mth temporal 
predictor on day j, and X 3m y is the value of the mth spatial-temporal 
predictor at site i on day j. gj(reg) and hj(reg) are the daily random 
intercepts and AOD slopes specific to each study area region, nested 
within the overall random effects Uj and vj. Here, we assume s is a 
3x3 diagonal matrix with diagonal elements S 2 U , S 2 V , S and s REG 
is a 2 x 2 diagonal matrix with diagonal elements S 2 g , S 2 h. 

We used ten-folds out-of-sample cross validation (CV) to vali- 
date our model 1 predictions. We randomly divide our data into 90 
and 10 percent splits ten times. We predict for the 10% datasets 
using the model fitted from the remaining 90% of the data. We then 
report these computed R 2 values. To test our results for bias we 
regress the measured PM value for a given site and day against the 
corresponding predicted value. We estimated the model prediction 
precision by taking the square root of the mean squared prediction 
errors (RMSPE). In addition we calculated prediction errors from a 
model that contained the spatial components only, to make it more 
comparable to the commonly used monthly/yearly prediction 
models available such as Yanosky et al. (2008). Temporal R 2 was 
calculated by regressing Delta PM against Delta predicted where: 
Delta PM is the difference between the actual PM in place i at time j 
and the annual mean PM at that location, and Delta predicted is 
defined similarly for the predicted values generated from the 
model. Spatial R 2 was calculated by regressing the site-specific 
annual means in observed PM versus the same annual means for 
predicted PM. 

The next model (model 2), uses the fit of model 1 to predict a 
PM 2.5 concentration for each day and grid cell at which we have an 
observed AOD value. This resulted in yearly datasets with PM 2.5 
prediction for all day-grid cell combinations with available AOD. 

In model 3 of the sequence, we estimate daily PM2.5 for all grid 
cells in the study area even in days when no AOD data are present. 
We fit a generalized additive model with a smooth function of 
latitude and longitude (using the grid cell centroids) and a random 
intercept for each cell. This is similar to other interpolation tech- 
niques (universal kriging, etc.) to use nearby grid cells to help fill in 
the missing, however we also construct a 100 km buffer around 
each grid cell and fit a regression analysis relating the predicted PM 
at each grid cell to the daily mean PM2.5 from the stations in each 
buffer around that grid cell, which also aids in predicting the value 
on missing days. The 100 km buffer size was chosen since we 
wanted a small enough buffer to ensure relevance and hence 
improve the R 2 of the monitored value on the grid cell values, but 
large enough a buffer to include multiple PM monitors to produce 
more stable estimates (again improving the prediction R 2 ). To allow 
for temporal variations in the spatial correlation, we fit a separate 
spatial surface for each two-month period of each year. Using this 
method provides additional information about the concentration in 
the missing grid cells that classic interpolation would not provide. 
Specifically, we fit the following semiparametric regression model: 

PredPMy = (a + u,) + (bi + v/jMPM,-/ + s(X,- ;Yj) fe ^ + y ; 

(W/v i) (00) ;U b ^ J 

where PredPMy is the predicted PM 2.5 concentration at a grid cell i 
on a dayj from the model 2; MPMy is the mean PM in the relevant 
100 km buffer for site i on a dayj; a and u,- are the fixed and grid-cell 
specific random intercepts, respectively: bj and v,- are the fixed and 


random slopes, respectively. The smooth X,-,Y,- are the latitude and 
longitude, respectively, of the centroid of grid cell i, and s(X,-,Y, ■)*<(/) is 
the a smooth function of location (modeled by thin plate splines) 
specific to the two-month period k(j) in which dayj falls (that is, a 
separate spatial smooth was fit for each two-month period). 

To estimate the goodness of fit and due to computational limits, 
we used a “leave 10 out” approach where we randomly selected 10 
monitors to leave out of the model 3 predictions. We the test the fit 
and bias (regressing the measured PM values against the predicted 
values and predicting PM levels at the left out monitoring 
locations). 

Finally, we run a forth model where we take the residuals 
constructed by taking the difference between a given daily moni- 
tored PM 2.5 concentration and the 1 km x 1 km corresponding daily 
model 3 prediction, and regress that against spatial and temporal 
predictors of PM 2.5 at each monitor. Specifically, we fit the following 
model: 

ResidPMy (Traffic density,- population density,-) 

+ f 2 ( Elevation,-) + / 3 (Percent urban,) 

+/ 4 (Distance to Al roads,-) 

+/ 5 (Distance to point emissions,-) +/ 6 PBL,-/ 

+/ 7 Traffic density, ;PBLy 

+/ 8 y Traffic density, Visibility y + y 

( 4 ) 

where ResidPMy is the residual at a spatial monitor site i on dayj]fi 
denotes a penalized spline for an interaction between traffic den- 
sity and population density at a spatial monitor site i; f 2 —f 6 denote 
(potentially nonlinear) effects of elevation, percent urban land use, 
distance to Al road, distance to point emissions and PBL respec- 
tively at a spatial monitor site i (and dayj for PBL). f 7 denotes a 
penalized spline for an interaction between traffic density and PBL 
at a spatial monitor site i on dayj and js denotes a penalized spline 
for an interaction between traffic density and visibility at a spatial 
monitor site i on dayj. Finally y is the error. 

3. Results 

Overall mean for MAIAC AOD was 0.25. Fig. 2 presents the mean 
ground PM 2.5 measurements across all years (2003—2011) and the 
difference between monitor and predicted PM 2.5 concentrations at 
each monitoring site. The figure shows how measured PM 2.5 con- 
centrations at monitors corresponds well with our predicted con- 
centrations, and the differences at 95% of the monitoring sites are 
within ±1.3 pg/m 3 , clearly showing excellent agreement between 
observed and predicted values. Fig. 3a presents a density plot 
exhibiting the daily variation of AOD slopes between 2003 and 2011 
during model 1 calibrations, while Fig. 3b presents a density plot 
exhibiting the daily variation of temperature slopes for the same 
period. The figure shows there is considerable day-to-day vari- 
ability in these slopes. 

Table 1 presents results from our model 1 calibration analysis. 
The yearly predictive power of our model is extremely high, with an 
overall “ out-of-sample ” R 2 for daily values of 0.88 (year-to-year 
variation 0.82—0.90), with a highly significant association between 
PM 2.5 and the main explanatory variable-AOD. The models yield 
almost no bias in our cross validation results as the slope of 
observed versus predicted = 0.99 (year-to-year variation 
0.98—1.01). The spatial and temporal out-of-sample results also 
presented excellent fits (Table 1): For the temporal model, the 
mean “ out-of-sample ” R 2 was 0.87 (year-to-year variation 



Fig. 2. EPA and IMPROVE PM 2 .5 annual means and the difference between measured and estimated PM2.5 concentrations at each PM2.5 monitor. 
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Daily AOD slopes Daily Temperature slopes 

Fig. 3. Density plots exhibiting the daily variation of AOD slopes (a) and temperature slopes (b) between 2003 and 2011 during the stage 1 calibrations. 


Table 1 

Prediction accuracy: Ten-fold cross-validated R 2 for PM 2 .5 stage 1 predictions (Calibration stage for 2003-2011). 


Year 

R 2 

Slope 

Spatial R 2 

Temporal R 2 

RMSPE (pg/m 3 ) 

Spatial RMSPE (pg/m 3 ) 

2003 

0.89 

0.98 

0.87 

0.88 

2.89 

1.05 

2004 

0.89 

0.99 

0.90 

0.88 

2.34 

0.85 

2005 

0.88 

0.98 

0.87 

0.87 

2.85 

1.00 

2006 

0.89 

1.00 

0.85 

0.88 

2.30 

0.81 

2007 

0.90 

1.01 

0.93 

0.90 

2.34 

0.73 

2008 

0.88 

0.99 

0.92 

0.86 

2.22 

0.72 

2009 

0.86 

0.98 

0.86 

0.84 

1.97 

0.71 

2010 

0.90 

0.99 

0.84 

0.89 

1.85 

0.69 

2011 

0.82 

0.99 

0.80 

0.81 

2.22 

0.82 

Overall Mean 2003- 

2011 0.88 

0.99 

0.87 

0.87 

2.33 

0.82 

0.81—0.90) and for the spatial model the mean “ out-of-sample ” R 2 

Fig. 5 shows the deviation of the estimated local PM 2.5 con- 

was 0.87 (year-to-year variation 0.80- 

0.93). When looking at Root 

centrations (model 4) from the average lxl km PM 2.5 concen- 

Mean Square Prediction Error (RMSPE), our models exhibit very 

trations at a very fine resolution (200 

x 200 m) aggregated over a 

low RMSPE values of 2.33 fig/nr (year-to-year 

variation 

year (2003) in Boston. 



1.95-2.89 |ig/m 3 ). 

Our “spatial” component only RMSPE is much 




lower at 0.82 pg/nr. The final prediction model (model 3) presents 

4. Discussion 



very similar results and is presented in Table 2. 





Fig. 4 shows the spatial pattern of predicted PM 2.5 

concentra- 

In this paper we used novel 1 km AOD MODIS data based on the 

tions in Boston from the AOD models, averaged over 

the entire 

MAIAC algorithm to 

predict PM 2.5 

concentrations across the 

study period. The spatial variation in 

these long term PM 2.5 pre- 

Northeastern USA. Using the newly available MAIAC data allows us 

dictions ranges from 2.36 pg/nr to 40.12 pg/m 3 showing a good 

to simplify our previously described models while still gaining 

range of variability for our model. 



much better predictive power and reducing the exposure error. In 


Table 2 

Prediction accuracy: R 2 for stage 3 PM 2.5 predictions (final prediction model including locations without AOD for 2003—2011). 

Year 

R 2 

Slope 

Spatial R 2 

Temporal R 2 

RMSPE (pg/m 3 ) 

Spatial RMSPE (pg/m 3 ) 

Slope for leave 10-out 

2003 

0.91 

1.02 

0.88 

0.91 

2.69 

0.65 

1.09 

2004 

0.89 

1.01 

0.93 

0.88 

2.47 

0.60 

1.09 

2005 

0.88 

1.01 

0.93 

0.88 

2.79 

0.63 

1.10 

2006 

0.89 

1.02 

0.93 

0.89 

2.37 

0.53 

1.08 

2007 

0.91 

1.01 

0.96 

0.91 

2.26 

0.47 

1.04 

2008 

0.87 

0.99 

0.96 

0.86 

2.30 

0.42 

1.03 

2009 

0.87 

1.01 

0.94 

0.86 

1.95 

0.37 

1.06 

2010 

0.89 

1.00 

0.96 

0.88 

1.91 

0.29 

1.04 

2011 

0.84 

1.01 

0.92 

0.83 

2.09 

0.39 

1.05 

Overall Mean 

0.88 

1.01 

0.93 

0.88 

2.32 

0.48 

1.07 
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Fig. 4. Mean PM 2 .5 concentrations in each 1 x 1 km grid cell in Boston during 2003 predicted by the AOD models. 


addition, in our updated model we introduce significant method- 
ological improvements with the predictions for daily fine resolu- 
tion (200 x 200 m) PM exposure, which can better capture 
pollution attributed to local sources such as traffic. Our models 
yield extremely good model fits and observed versus predicted 
slopes that were practically T, which clearly show that there is no 
bias in the resulting exposure predictions. 

It is important to compare the predictions of our models to the 
other ‘state of the art’ prediction models that are available today. 


While the various existing models have all advanced in the past few 
years (including some models that use MAIAC data as well), they all 
still lack the ability to generate daily predictions for studies of the 
acute effects (short term) of air pollution as well as chronic (long 
term) effects, which are useful in epidemiologic studies aiming to 
estimate both acute and chronic effects of exposure. Sampson et al. 
(2013) recently published a universal kriging model using Partial 
Least Squares regression for estimating annual PM2.5 concentra- 
tions. They present a high accuracy of prediction with an overall R 2 



Fig. 5. The deviations of the estimated local pollution concentrations at a very fine resolution (200 x 200 m) from the average 1 x 1 km grid PM2.5 aggregated over a year (2003) in 
Boston. 
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of 0.88 and well-calibrated predictive intervals. While the predic- 
tion accuracy of this model for long-term exposures is high, this 
model ( 1 ) does not predict daily and ( 2 ) is complex and requires an 
enormous number of spatial variables that are not always freely- 
available to researchers. In contrast our model presents similar 
spatial accuracy ( R 2 = 0.87), but with the additional advantage of 
yielding daily predictions also with very good prediction accuracy 
(temporal R 2 = 0.87). Hu et al. (2014) published a study estimating 
ground-level PM 2.5 concentrations in the Southeastern United 
States. Similar to our model they used MODIS based MAIAC data 
with a two-stage spatial statistical model with meteorological 
fields and land use parameters as ancillary variables to estimate 
daily mean PM 2.5 concentrations. The model was run in south- 
eastern USA, for the year 2003. They used cross validation and 
predicted for only days with available MAIAC data resulting in an R 2 
of 0.67 and RMSPE of 3.88 pg/m 3 . Their models do not yield daily 
predictions (in every grid cell on every day) and are less accurate 
compared to our presented model (such as cross validated RMSPE 
of 3.88 pg/m 3 compared to 2.33 pg/m 3 in our model). Other models 
published recently have all similar issues with lack of highly ac- 
curate daily predictions across the study area (Kim et al., 2013; Lin 
et al., 2013; Liu et al., 2009; Yanosky et al., 2008; Yap and Hashim, 
2013). It is important to note that Beckerman (Beckerman et al., 
2013) recently published results showing that in a Hybrid model 
estimating national scale spatio-temporal variability of pni 2.5 the 
inclusion of satellite based AOD into the LUR model did not improve 
the model. The study though used the older 10 x 10 MODIS data 
and generally the cross validated results were lower than our 
presented models (R 2 = 0.79). 

Although the new lxl km MAIAC data is a considerable 
improvement over the previous MODIS 10 x 10 km data (for 
example, R 2 = 0.83 for our previously published 10 km model 
versus R 2 = 0.88 in the current MAIAC model, spatial R 2 = 0.78 
versus R 2 = 0.87), finer scale or less noisy future satellite-based data 
could further reduce exposure error bias resulting in larger and 
more accurate health effects estimates. Jerrett et al. (2005) have 
previously demonstrated how fine-scale variations in PM 2.5 are 
associated with larger health effects than those that vary regionally. 
We partly address this limitation in our model by adding our new 
local stage model 4 predictions at very fine 200 x 200 m scale that 
complements our lxl km predictions. 
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